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ABSTRACT 

It is widely accepted that cosmic rays (CRs) up to at least PeV energies are Galactic in 
origin. Accelerated particles are injected into the interstellar medium where they propagate to 
the farthest reaches of the Milky Way, including a surrounding halo. The composition of CRs 
coming to the solar system can be measured directly and has been used to infer the details of CR 
propagation that are extrapolated to the whole Galaxy. In contrast, indirect methods, such as 
observations of 7 -ray emission from CR interactions with interstellar gas, have been employed to 
directly probe the CR densities in distant locations throughout the Galactic plane. In this article 
we use 73 months of data from the Fermi Large Area Telescope in the energy range between 
300 MeV and 10 GeV to search for 7 -ray emission produced by CR interactions in several high- 
and intermediate-velocity clouds located at up to ^ 7 kpc above the Galactic plane. We achieve 
the first detection of intermediate-velocity clouds in 7 rays and set upper limits on the emission 
from the remaining targets, thereby tracing the distribution of CR nuclei in the halo for the first 
time. We find that the 7 -ray emissivity per H atom decreases with increasing distance from the 
plane at 97.5% confidence level. This corroborates the notion that CRs at the relevant energies 
originate in the Galactic disk. The emissivity of the upper intermediate-velocity Arch hints at a 
50% decline of CR densities within 2 kpc from the plane. We compare our results to predictions 
of CR propagation models. 
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1. INTRODUCTION 

Cosmic rays (CRs) are a wide-spread phe¬ 
nomenon, pervading galaxies and inter galactic 
space. They are directly detected within the so¬ 
lar system with spacecraft, balloon-borne, and 
ground-based instruments. Other methods to 
study CRs include observations of radio syn¬ 
chrotron emission and 7 rays produced in interac¬ 
tions of CRs with interstellar gas, radiation and 
magnetic fields. CRs consist of all known sta¬ 
ble and long-lived isotopes ranging from the most 
abundant hydrogen and helium nuclei through the 
traces of very rare Actinides, and particles such as 
electrons, positrons, and antiprotons. The energy 
density of CRs in the interstellar medium (ISM) 
is comparable to those of the Galactic interstellar 
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radiation field, magnetic field, and turbulent mo¬ 
tions of the interstellar gas. CRs are one of the es¬ 
sential components of the Galactic ecosystem (e.g., 
Ferriere|2QQl ), affecting the thermal, chemical and 
magnetohydrodynamical state of the ISM (e.g., 
Ptuskin et al. 2QQ6[ [Indriolo et al.|[2QQ^ . Emis¬ 
sion from CRs interacting within the Milky Way 
is also a source of foregrounds/backgrounds for 
many observations spanning from radio to 7 rays. 

Observations of X-ray and 7 -ray emission from 
Galactic objects such as supernova remnants, pul¬ 
sars and massive-star clusters reveal the presence 
of energetic particles, suggesting that efficient ac¬ 
celeration processes occur in their vicinities (for re¬ 
cent reviews see, e.g., Reynolds 2008 Bykov|2014 ). 
The Galactic origin of CRs is also supported by 7 - 
ray observations of external galaxies, which indi¬ 
cate that the CR densities are not universal (e.g., 
Sreekumar et al. 1993) and are correlated with 


the star-formation properties of the galaxies (e.g., 
Ackermann et al^| 2012 b ). Also, the radial sizes of 
radio halos measured in external, edge-on spiral 
galaxies appear to be correlated with active star 
formation in their disks ( Dahlem et al.|p!9^ ). 

CR particles accelerated in sources in the 
Galactic disk are injected into the ISM where 
they propagate through a surrounding halo before 


escaping into the intergalactic space (e.g.. Strong 


et al. 2010). This paradigm is substantiated in 


models that can reproduce reasonably well all the 
related observables (for a recent review see, e.g.. 
Strong et al. 2007[ ). The propagation of particles in 
the ISM leads to the destruction of primary nuclei 
via spallation and gives rise to secondary nuclei 
and isotopes which are rare in nature, and other 
particles. The observed abundances of stable sec¬ 
ondary CR nuclei (e.g., boron) and radioactive iso¬ 
topes (^^Be, ^^Al, ^^Cl, ^^Mn) allow the separate 
determination of the confinement halo size and 
diffusion coefficient (e.g., [Strong fc Moskalenko 


1998) which is somewhat model dependent. Typ¬ 


ical constraints on the halo height arising from 


such method are 4-6 kpc (e.g., Moskalenko et al. 
^otta et ar]|2Qll ). 
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Observations of edge-on external galaxies show 
directly the existence of radio halos extending 
to kpc distances perpendicularly to their disks 


(e.g., Dahlem et al. 1994), and sometimes up to 
^ 15 kpc ( Irwin et al.||1999 ). Radio observations 
of large-scale synchrotron emission from the Milky 


Way were also interpreted in the context of diffu¬ 
sion models as indicating a vertical extent of the 
halo of 4 kpc to 10 kpc (e.g.. Strong et al.||2QQQ ), 
with values of ^ 10 kpc favored based on a re¬ 
cent analysis including WMAP and 408 MHz ob¬ 
servations ( Orlando fc Strong||2013 ). The discrep¬ 
ancies between halo sizes determined with differ¬ 
ent methods are not necessarily surprising. The 
method based on the direct detection of GR nu¬ 
clei provides the size of the effective volume filled 
with GRs that are still confined to the Milky Way. 
In contrast, observations of synchrotron emission 
show larger halos because part of the emission may 
be generated by particles that are leaking into the 
intergalactic space, i.e., that are already beyond 
the “point of no return”. 

The Large Area Telescope (LAT) on the Fermi 
Gamma-ray Space Telescope mission (Atwood 
et al.||2QQ9 ) has been surveying the sky since 2008 


in a wide energy range of 20 MeV to > 300 GeV. 
These observations provide information about GR 
particles (mostly protons and electrons) in distant 
locations that can be used to study CR propaga¬ 
tion (e.g., [Ackermann et al.| 2012 a| ). Measurements 
with the LAT confirmed earlier observations (e.g.. 
Strong et al.[[198^ that the radial gradient of the 


7 -ray emissivity ( 7 -ray emission rate per hydro¬ 
gen atom) in the ISM outside the solar circle is 
much smaller than that in the density of putative 
CR sources ( [Abdo et al. 2010 Ackermann et al 


2011). This can be interpreted as another hint of 


a ^ 10 kpc CR confinement halo, because it facil¬ 
itates radial diffusion of CRs and compensates for 
the paucity of CR sources in the outer Galaxy. 

In this work we propose a more direct method 
to constrain the propagation of CRs in the halo 
of the Milky Way with LAT data. Although the 
interstellar gas in the Milky Way is largely con¬ 
fined to the Galactic plane with a scale height for 
atomic hydrogen of < 300 pc across much of the 
disk ( Kalberla fc Ker^[2009 ), certain populations 
of interstellar clouds are known to exist at much 
greater distances from the plane. So-called high- 
velocity clouds (HVCs) and intermediate-velocity 
clouds (IVCs) are interstellar clouds of atomic hy¬ 
drogen, Hi, with line-of-sight velocities inconsis¬ 


tent with Galactic rotation (see, e.g., van Woerden 
et al.|20Q5 ). The demarcation between HVCs and 


IVCs is conventionally taken to be Doppler shifted 
velocity with absolute value between 80 km s“^ 
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and 100 km s ^ with respect to the local standard 
of rest. 


The distances and metallicities of HVCs/IVCs 
can be inferred from spectroscopic observations 
of foreground and background stars. HVCs typ¬ 
ically have metallicities less than Solar and are 
generally considered to be remnants of the for¬ 
mation of the Milky Way, having existed in the 
halo of the Galaxy since the disk formed, or old 
tidal streams extracted from nearby dwarf galax¬ 
ies (Wakker 2002). Some IVCs have metallici¬ 
ties much closer to Solar and these may be gas 
ejected from the plane, e.g., by Galactic fountains 
( Wakker||2002 ). No molecular gas, e.g., CO emis¬ 


sion, has been detected in HVCs (Wakker et al. 


1997), and IVCs are also typically free of CO emis¬ 


sion as well. 


These clouds have low column densities, V(H i) ^ 
10i9 — 10^^ cm“^, but can subtend hundreds of 
square degrees. The exposure of the LAT sky sur¬ 
vey is deep enough for HVCs/IVCs to conceivably 
be detected in diffuse y-ray emission and used to 
estimate the CR densities in the Galactic halo. 
These clouds are not actively forming stars and 
are believed to be free of internal sources of CRs 
and hence any y-ray emission from them samples 
the large-scale distribution of Galactic CRs. 

The article is structured as follows. In Section[2] 
we describe our criteria for selecting a sample of 
HVCs/IVCs for study. In Section we discuss 
the y-ray and multi wavelength survey data we use. 
Section [^provides the details of how we derive the 
Hi maps, in particular how we fit spectral com¬ 
ponents with different line-of-sight velocities to 
achieve kinematic separation of interstellar gas by 
distance along the line of sight, and other maps of 
ISM tracers to accurately account for foreground 
gas in the local neighborhood. In Section we in¬ 
troduce the y-ray analysis procedure to measure 
or constrain the CR densities in the targets, and 
the procedure we implemented to estimate sys¬ 
tematic uncertainties. In Section |6] we discuss the 
results for the individual regions that we studied 
and the implications for the propagation of CRs in 
the halo of the Milky Way. The conclusions are in 
Section Appendix [^provides more technical in¬ 
formation about the derivation of the maps of the 
ISM, including an iterative approach to unbiasing 
the maps of the column density in the dark neu¬ 
tral medium (DNM). Appendix [B] describes how 


we computed uncertainties and upper limits with 
the jackknife method for evaluating the systematic 
uncertainties. 


TARGETS 


Wakker (2001) reports the most complete cen¬ 


sus to date of IVCs and HVCs, including lower 
and upper limits on their distances (which we refer 
to as distance brackets), from spectrophotometric 
distances of foreground and background stars. In 
order to be able to constrain the CR densities in 
well-defined regions of the Milky Way halo, we se¬ 


lect as targets for our analysis clouds in [Wakker 
(2001 | ) for which distance (altitude) brackets have 
been determined. We summarize the selected tar¬ 
gets in Tablegrouped in three regions of interest 
(ROIs) as used for the analysis later. In each case 
the ROIs enclose the parts of the clouds associated 
with the foreground and background stars used for 
determining the distance brackets. 

Region A encompasses an IVC known as the 
low-latitude intermediate velocity (IV) Arch, and 
an HVC known as Complex A. The low-latitude 
IV Arch has a mass of 1.5 — 6 x 10^ Mq at a dis¬ 
tance of 0.9 — 1.8 kpc from the Sun, corresponding 
to a distance above the plane of 2 ; = 0.6 — 1.2 kpc. 
According to Wakker | ( |2001 ) the distance (compa¬ 
rable to that of the Perseus arm of the Galaxy), 
the metallicity (close to solar), and negative veloc¬ 
ities in the range of —20 km s“^ to —30 km s“^ 
(higher than the Perseus arm itself in this direc¬ 
tion) are all compatible with the characteristics 
of a high-interarm cloud that is part of the re¬ 
turn flow of gas injected from the disk into the 
halo by a fountain process. Complex A is a dis¬ 
tant HVC with the narrowest distance brackelElof 
4.0 — 9.9 kpc {z = 2.6 — 6.8 kpc), which implies 
a mass of (0.3 — 2) x 10^ Mq. Its low metallic¬ 
ity suggests that this HVC is a relic of the ancient 
Universe, a clump of gas that never underwent star 
formation and is now falling into the Milky Way. 

Region B encompasses two IVCs that are 
together known as the IV Arch. The gas in 
the IV Arch is seen in two distinct velocity 
ranges: > —60 km s“^ (the lower IV Arch) and 


^HVC Co mplex C has a broad er distance bracket of 3.7 — 
11.2 kpc ( [Wakker et al.|2007| . We leave it for future stud¬ 
ies. 
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Table 1 

Target regions and gomplexes 


Region 

ROI 


Complex 

d 

z 

Metallicity 

^^LSR 


Center (Z, b) 

Pixels 


(kpc) 

(kpc) 

(relative to solar) 

(km s^L 

A 

(155°,37°) 

161,81 

Low-Latitude IV Arch 

0.9-1.8 

0.6-1.2 

1.0 ± 0.5 

—90 < vlsr < —40 




Complex A 

4.0-9.9 

2.6-6.8 

0.02-0.4 

V-LSK < —90 

B 

(182°,57°) 

125,141 

Lower IV Arch 

0.4-1.9 

0.4-1.7 

~ 1 

—60 < vlsr < —15 




Upper IV Arch 

0.8-1.8 

0.7-1.7 

~ 1 

^^LSR < -60 

C 

(251°,69°) 

85,117 

IV Spur 

0.3-2.1 

0.3-2.1 


—90 < vlsr < —20 


Note. —The ROIs are defined as rectangular regions in a plate carree projection in Galactic coordinates, centered at each ROI 
center position and with the number of pixels of 0?25 width given above. See text for details. 

Note. —The brackets in d and 2;, and the metallicities are taken from [Wakker||200l| ) . The metallicity of the IV Spur is unknown. 

Note . —T he ranges in velocities with respect to the local standard of rest, vlsr, are the preliminary boundaries that we use in 
Section [4.1| to construct the V(Hi) maps. 


< —60 km s“^ (the upper IV ArclQ. We consider 
in this analysis only the portion of the IV Arch 
seen at Galactic longitudes I > 150°, for which 
there is a distance bracket for the lower part of the 
Arch of 0.4 — 1.9 kpc {z = 0.4 — 1.7 kpc). At lower 
longitudes the distance is less reliably determined, 
and part of the lower IV Arch may be associated 
with a CO core with distance < 0.3 kpc from 
the solar system ( [Wakker 2001[ ). On the other 
hand, the upper IV Arch has a distance bracket 
of 0.8 — 1.8 kpc {z = 0.7 — 1.7 kpc). Therefore, 
whether the two velocity components of the IV 
Arch are also physically separated entities is un¬ 
certain. The total mass of the IV Arch is of the 
order of several 10 ^ Mq, and its metallicity is 
not well determined, but possibly close to solar 
( |Wakker|[2M ). 


Region C encompasses the IV Spur, which is an 
extension of the IV Arch at lower velocities and 
higher Galactic latitudes. The distance bracket is 
0.3 — 2.1 kpc {z = 0.3 — 2.1 kpc), which implies 
a mass of 0.2 — 2.8 x 10^ Mq. The metallicity is 
unknown. 


All the target objects discussed reside in a rel¬ 
atively small sector of the Milky Way with 130° < 
I < 280° and b > 25°. This is in fact one of the 
in IVGs and HVGs, and for which 
could get access to stellar probes 

observations. 


regions richest 


Wakker (2001 


^The nomenclature refers to the absolute value of the veloc¬ 
ity. 


Furthermore, in recent years distance brackets 
have been determined for additional, lower-mass 
IVGs and HVGs (e.g., Wakker et al.|[2QQ8 ). Such 


clouds have ratios of mass over distance square 
of the order of < 10^ Mq kpc“^. Therefore, for 
normal y-ray emission rates per H atom, as mea¬ 
sured in the Galactic disk, we expect total fluxes of 


< 10“^^ cm“^ s“^ at energies > 100 MeV. These 
fluxes are one order of magnitude or more below 
the LAT sensitivity level (Ackermann et al. 2012c). 
Therefore, we do not expect to be able to use LAT 
observations of these clouds to constrain their CR 
densities, and we do not include them among our 
targets. Gonversely, the clouds selected for the 
analysis all have mass over distance square ratios 
of the order of > 10 ^ Mq kpc“^, and hence ex¬ 
pected 7 -ray fluxes > 10“^ cm“^ s“^, in the range 
detectable by the LAT. 

We have not included among our targets Gom- 
plex GGP, also known as the Smith Gloud, because 
this HVG is seen at Galactic latitudes \b\ < 20° 
toward the inner Galaxy, where the diffuse y-ray 
foregrounds from the Galactic disk are brighter 
and more challenging to model. We leave Gom- 
plex GGP for future studies. We note that |Drli^ 


Wagner et al. (2014) and Nichols et al. (2014) 


searched for y-ray emission associated with Gom- 
plex GGP due to exotic processes related to hypo¬ 
thetical dark matter particles, and did not report 
on any significant detections. 
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3. DATA 


3.1. 7 - Ray Data 

The LAT is a pair-conversion 7 -ray telescope 
for the energy range 20 MeV to > 300 GeV. The 


instrument is described in Atwood et ah (2009), 
and its on-orbit performance discussed in [Acker- 1 
mann et ah (2012c). We use LAT observations 


accumulated from the beginning of the Science 
phase of the mission, METQ 239500801 (2008 Au¬ 
gust 4), up to MET 431481603 (2014 September 
4), for a total of 73 months. We use the P7REP 
dataset ( jBregeon et al. 2013) and select “Clean” 
class events in order to limit the contamination 
from backgrounds due to residual CRs misclassi- 
fied as 7 rays in the study of large-scale diffuse 
features. 

We apply further restrictions to the dataset 
used: 


• We select only time intervals when the LAT 
was operated in normal Science mode (e.g., 
excluding calibration data) and when the 
data quality was considered good (e.g., ex¬ 
cluding solar flares which caused pile-up phe¬ 
nomena in the detector). 

• We retain only candidate 7 rays detected at 
inclinations from the Earth zenith < 95°, in 
order to limit the contamination from the 
bright 7 -ray emission produced by CR inter¬ 
actions in the Earth’s atmosphere that forms 
extended features in the sky. 


We consider only candidate 7 rays with 
measured energies between 300 MeV and 
10 GeV. Below 300 MeV the LAT angular 


resolution rapidly deteriorates (Ackermann 


et al. 2012 c), which makes separating dif¬ 


ferent components of diffuse 7 -ray emission 
based on their morphologies more difficult. 
Additionally, below ^ 300 MeV one would 
need to model the residual emission from CR 
interactions in the upper atmosphere. This 
can introduce a faint diffuse glow at moder¬ 
ately high declinations, which would need to 
be modeled, and would further complicate 
the study of faint extended features espe¬ 
cially near the North Celestial pole where 


^ Fermi Mission Elapsed Time, measured in seconds since 1 
January 2000. 


our targets are located. Above 10 GeV the 
photon counts are sparse and the data are 
not very useful to constrain CR interactions 
in clouds owing to their spectrum being 
softer than inverse-Compton (IC) emission 
from Galactic CR electrons and isotropic 
backgrounds (e.g., Ackermann et al.||2Q12a ). 
Limiting the energy range also helps to at¬ 
tenuate the systematic uncertainties due to 
the spectral models assumed for the various 
components of the 7 -ray sky. 

Eor the 7 -ray analysis described in this section 
and elsewhere in the article we used the official 
LAT Science Tools^ version v09r34p0l[^ 

3.2. Interstellar Medium Tracer Data 


If the CR densities are uniform within a cloud, 
the 7 -ray intensities are simply proportional to the 
target gas column densities. Thus, in order to 
model 7 -ray emission from IVCs and HVCs, and 
foreground emission from local gas we rely on a 
set of ISM tracers described here. 

We employ as tracer of atomic gas the Leiden- 
Argentine-Bonn (LAB) all-sky survey of the 21 - 
Hi emission line by [Kalberla et al 


cm ELI emission line by [Kalberla et al.| (2005). 
The survey provides the Hi line brightness tem¬ 
peratures over a grid of 0?5 (with an effective an¬ 
gular resolution of ^ 0?6), and covers velocities 
with respect to the local standard of rest from 
—450 km s“^ to +400 km s“^, at a resolution 
of 1.3 km s“^. Data were corrected for stray ra¬ 
diation reaching a residual noise contamination 
< 0.09 K for most directions in the sky. We use 


data from Kalberla et al. (2005) to derive V(Hi) 
column-density maps as described in Section [TTj 
Molecular hydrogen cannot be traced directly 
in its most common cold state. The most widely 
used surrogate tracer is the 2.6 mm line of ^^CO. 
Eor this work we use a map of CO brightness 
temperature integrated over Doppler velocities, 
Wco^ derived from the Center for Astrophysics 
composite CO survey (Dame et al. 2001 ) supple¬ 
mented with higher-latitude observations of the 


Ursa Major region (de Vries et al.[1987). The CO 
data, denoised using the moment-masking tech¬ 


nique (Dame 2011), were used to construct the 


"^The Science Tools are available along with LAT data from 
the Fermi Science Support Center, http://fermi.gsfc. 
nasa.gov/ssc 
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Wco iiiap on a 0?25 grid. The VFco intensities 
are assumed to be proportional to the molecular 
hydrogen column densities A/'(H 2 ). 

Only ROI A contains significant CO emission 
(Figure [^. The CO emission detected is as- 
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Wolfire et al.||2QlQ Smith et al.||2Q14 ). There can 
also be a contribution from atomic gas overlooked 
because of the approximations applied in deriv¬ 
ing column densities from the Hi line intensities 
(Fukui et al. 2015). On large scales the DNM 


in the Milky Way is best traced by dust ther¬ 
mal emission or extinction and 7 rays, which re¬ 
veal correlated excesses on top of the components 
traced by Hi and CO that can be interpreted as a 
contribution from undetected, but otherwise nor¬ 


mal, neutral interstellar matter (e.g., Grenier et al 


2QQ5||Abdo et al.|2QlQ[ |P/ancA; and Fermi Collab- 


orations ||2014| | 

In order to trace the DNM in this work we 
rely on the all-sky model of thermal dust emission 
by the |P/ancA; Collaboration XI | (2Q14) based on 
Planck and IRAS data {Planck public data release 
R 1 . 20 ). We test both the map of optical depth at 
353 GHz, r 353 , and the map of total dust radi¬ 
ance, R. The angular resolution of these maps is 
5'. They are used to construct maps of the DNM 
as described in [4^ 


Fig. 1 .— Wco intensities in the direction of 
ROI A. The map is shown in zenithal equidistant 
projection. Blank pixels are set to zero. The black 
line shows the border of the ROI. 


sociated with the Ursa Major molecular clouds, 
located at a few hundred parsecs from the so¬ 
lar system (e.g., de Vries et al. 1987). Note 


that, according to the results from the component 
separation of the all-sky survey performed using 


the Planck-UFl instrument (Planck Gollaboration 
XIII I2014D , there is no significant GO emission in 
our ROIs besides what is included in this CO map. 

Various observations indicate that a linear com¬ 
bination of V(Hi) and Wco maps does not prop¬ 
erly account for the totality of the ISM neutral 
gas, but large fractions of the gas can be in a 
DNM phase, also known as dark gas (e.g., Mag- 


et al. 


2QQ3[ Grenier et al. 2QQ5[ |Abdo et~ 


2QIQ| Langer et al. 2QIQ| Planck and Fermi Gol 


laborations 2014). The DNM may be composed 


of molecular gas not associated with GO, revealed. 


e.g., by alternative tracers like CH (e.g., Magnani 
et al.||2QQ3 ), that is likely to exist because the CO 
molecule is more prone to photodissociation than 
H 2 in the external layers of molecular clouds (e.g.. 


4. CONSTRUCTION OF THE INTER¬ 
STELLAR MEDIUM TRACER MAPS 

For each ROI we construct maps of the V(Hi) 
column densities for different complexes along the 
line of sight, and maps of the DNM column densi¬ 
ties. The maps are constructed over regions that 
include at least 5° boundaries around the ROI 
in order to properly take into account the point- 
spread function (PSF) of the LAT in fits of the 
models to the 7 -ray data. 


4.1. Construction of the Hi Maps 


The column densities of atomic gas, V(Hi), in 
IVCs and HVCs, and those in the local foreground 
gas can be separated based on Doppler velocities in 
order to separate contributions from CR interac¬ 
tions with gas in different locations in the Galaxy. 
To do so we adapt the procedure described in 


Abdo et al. (2010). The procedure used in this 


work has three steps. 


i For each ROI we define some preliminary 
velocity boundaries that separate differ¬ 
ent components along the line of sight, 
namely low-velocity (local) gas, and gas in 
HVCs/IVCs (Table g. 
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ii For each line of sight on the grid of the LAB 
survey we adjust the velocity boundaries so 
that they fall at the closest minimum (or 
inflection point) in the brightness temper¬ 
ature profile. The column density A/'(Hi) 
associated with each component is propor¬ 
tional to the integral of the brightness tem¬ 
perature over velocity within the adjusted 
velocity boundaries. 

iii For each line of sight on the Hi grid we fit 
to the brightness temperature profile a com¬ 
bination of Gaussian functions centered at 
the profile maxima (or inflection points if a 
maximum is not found for a velocity compo¬ 
nent) and use the results of the fit to correct 
the column densities of each component for 
the spillover from/into the adjacent compo¬ 
nents. 


The procedure is illustrated for an example line of 
sight in Figure Technical details of the proce¬ 
dure are illustrated by Tibaldo (2011 2.1.2, Ap¬ 
pendix B). 

In region A we set the preliminary velocity 
boundaries at —40 km s~^ to separate local gas 
from the low-latitude IV Arch, and at —90 km s“^ 
to separate the low-latitude IV Arch from Com¬ 
plex A. At latitudes b < 24° there is a significant 
contribution at higher velocities from gas in the 
disk of the Galaxy belonging to the Perseus and 
outer arms. Because Abdo et al. (2010) found that 
the 7 -ray emissivity of the gas in these arms is 
compatible with the local emissivity, we assign all 
of this gas to the local component. We note that 
in any case this low-latitude region is not included 
in the ROI, but only used as part of the border to 
ensure a proper convolution of the model with the 
LAT point-spread function. 


In region B we set the preliminary velocity 
boundaries at —15 km s~^ to separate local gas 
from the lower IV Arch, and at —60 km s“^ to 
separate the lower from upper IV Arch. In the 
case of region C we set the preliminary velocity 
boundaries at —20 km s“^ to separate local gas 
from the IV Spur, and at —90 km s~^ to separate 
the IV Spur from high-velocity gas; no significant 
emission is present at high velocities and the cor¬ 
responding map is not used in the construction of 
the 7 -ray model. 




VsRlknns M 


Fig. 2.— Hi brightness temperature, Tb, as a 
function of Doppler velocity with respect to the 
local standard of rest, Vlsr, for the line of sight at 
I = 172?5, b = 70?5. The upper panel illustrates 
how the preliminary velocity boundaries (dotted 
vertical lines) are adjusted to fall at the near¬ 
est minima or inflection points of the Tb profile 
(dashed vertical lines). The bottom panel illus¬ 
trates the fitting procedure to correct the V(Hi) 
column densities for spillover from one region to 
the next: in red Gaussian functions associated 
with peaks of the IV component, in green the 
Gaussian function associated with the peak of the 
low-velocity (local) component, in blue the total 
(no peaks were detected in this direction for the 
high-velocity component). 
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The Gaussian fits reproduce the measured 
brightness temperature profiles within ^ 5% on 
average (for the ratio of the integral of the ab¬ 
solute values of the differences between measured 
and fitted temperature along the line of sight to 
the integral of measured temperatures). The fit 
is significantly worse along a few line^ of sight 
with differences up to ^ 25%. The magnitude of 
the correction for the spillover is on average 3% of 
the uncorrected column densities. (The average is 
the largest for ROI B, where it is ^ 6%.) Even 
though the correction is only approximate, ow¬ 
ing to the imperfect modeling of the temperature 
profiles by the Gaussian functions, the average 
magnitudes of the errors are < 0.5% of the total 
column densities. In a small subset of the lines 
of sight (^ 0.5%) the fits do not converge prop¬ 
erly, mostly due to the faintness of H i emission in 
these high-latitude regions. In such cases we do 
not apply any correction for spillover. 

We have calculated the A^(H i) column densities 
under the approximation of small optical depth. 
An approximate correction for Hi optical depth 
based on the assumption of a uniform spin tem¬ 
perature has often been used in the literature (e.g.. 



regions studied make the effect of the correction 
small. Assuming a uniform spin temperature of 


80 K (as in 

would increase the column densities on average at 
most by 5% for the local component in region A, 
and less for other components and regions. We 
also note that, by using dust maps to trace gas 
missed by the A/'(Hi) and Wqo maps as described 
in |4.2[ we are effectively accounting also for H i 
densities underestimated due to the approxima¬ 
tion of small optical depth. 

We replace the N(Ri) column density seen in 
ROI A at / = 138.5°, b = 37° with the average of 
neighboring lines of sight owing to an anomalous 
temperature profile indicative of an artifact in the 
LAB data. 

The resulting N(Ri) maps are shown in Fig¬ 
ures i B and Figure (a) and (b) show 
some level of anticorrelation between higher/lower 
A/'(Hi) intensity regions. This is due to some gen- 


Planck Gollaboration XXIV 2011 


^Out of 11421, 34001, 32421 analyzed in region A, B, and 
C, respectively. 


uine anticorrelation between the intensities of H i 
emission in velocity space. At the same time, 
the lower IV Arch shows a bright feature at ^ 
—30 km s“^ that blends into the local Hi com¬ 
ponents, which makes the separation particularly 
challenging. This is addressed by the spillover cor¬ 
rection. 


4.2. Construction of the DNM Maps 


Since the work of Grenier et al. (2005) on the 


construction of models of interstellar y-ray emis¬ 
sion, dust residuals, i.e., total dust opacity or 
extinction minus the best-fit linear combination 
of Hi and Wco maps, have been used to trace 
the DNM, neutral interstellar matter not properly 
traced by the linear combination of V(Hi) and 
Hco- However, the original implementation of 
this method has two limitations: 1) the fit of the 
dust map is biased by the missing DNM compo¬ 
nent; 2) the effect of the assumptions on the stated 
errors in the fit of the dust maps is overlooked. 
Recently, th.e \Planck and Fermi Gollaborations 


[ (2014) performed an analysis of Planck and LAT 
data of the Chameleon interstellar complex where 
they refined the original method by: 1) developing 
a “circular” fitting method, where dust residuals 
are used to model the DNM in the y-ray analy¬ 
sis, and y-ray residuals to model the DNM in the 
dust analysis iteratively until a stable solution is 
found; 2) exploring several assumptions on the er¬ 
rors used in the fits and folding the differences into 
the errors on the final results. 

The method developed by the \Planck and| 


Fermi Collaborations (2014) derives a map of 


the DNM column densities from y rays to use in 
fitting the dust maps. Owing to the low gas col¬ 
umn densities (on average one order of magnitude 
lower than in the Chameleon complex) and the 
limited statistics of the y-ray dataset, it was not 
applicable to our ROIs. We therefore developed 
an alternative method to acknowledge and address 
the issue of the missing DNM component in the 
dust fit that does not rely on y-ray data and is 
described in the rest of the section. 

Owing to the limited angular resolution of y-ray 
data and the other ISM tracers, and to the input 
projections accepted by the LAT Science Tools^ 


we reproject the maps from the Planck Collabo¬ 


ration XI (2014) onto the grid chosen for the y-ray 
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Fig. 3.— A^(Hi) column density maps for ROI A: a) Low-velocity (local) component, b) the low-latitude 
IV Arch, c) Complex A. The black lines show the border of ROI A, as used in the analysis. The maps are 
shown in zenithal equidistant projection. 
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Fig. 4.— A^(Hi) column density maps for ROI B: a) Low-velocity (local) component, b) the lower IV Arch, 
c) the upper IV Arch. The black lines show the border of ROI B, as used in the analysis. The maps are 
shown in zenithal equidistant projection. 
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Fig. 5.— A/'(Hi) column density maps for ROI C: a) Low-velocity (local) component, b) the IV Spur. There 
is no significant Hi emission at high velocities in the ROI. The black lines show the border of ROI C, as 
used in the analysis. The maps are shown in zenithal equidistant projection. 


analysis in each ROI. 

Under the assumptions that dust grains are well 
mixed with gas in the cold and warm phases of the 
ISM (e.g., Bohlin et al.||1978 ), and that the prop¬ 
erties of dust grains and photon fields are uniform 
within each phase and complex, we build a model 
for the dust map D (i.e., the map of optical depth 
at 353 GHz, rsss, or radiance, R) by assuming that 
either quantity is proportional to the atomic and 
molecular column densities as traced by V(Hi) 
and Wco^ respectively, plus an isotropic term rep¬ 
resenting any zero-level shift, for example, due to 
residual contamination from the cosmic infrared 
background. 


3 

M°{l,b) = Y,ymymi)il,b)i+y°co-Wco{l,b)+Di, 
1=1 

( 1 ) 

Here and in the following equations we indicate 
with the dust specific power (opacity) per 

hydrogen atom, and with yco the specific power 
(opacity) per Wqo unit in the molecular phase. 

We use indices i = 1,2,3 to label the complexes 
separated kinematic ally along the line of sight as 
described in Section ItTI The model is fit to 
the dust map D using the minimum method. 


The is defined, for a generic model M, as 


l,b 


[Dii,b) - Mii,b)y^ 
a^(l,b) 


( 2 ) 


Following the work of the Planck and Fermi Col- 
laborations (2014), we consider different possi¬ 


ble assumptions on the uncertainties a{l^b). Our 
baseline assumption is that the uncertainties are 
mainly due to imperfections of the model because 
of the assumption of proportionality between ra¬ 
diance/opacity and gas column densities. In the 
absence of a reliable estimate of the uncertain¬ 
ties on M, we assume the uncertainties to be 
cr(/,5) oc M(/, 6 ). We empirically consider two al¬ 
ternative hypotheses to assess the impact of this 
choice: 


• cr(/, 6) (X !)(/, 6), which tests how the agree¬ 
ment between M and D affects the results 
of the fits; 

• cr(/,6) oc where (Jr{l^h) is the map 

of statistical uncertainties on rsss provided 
by the |P/ancA; Collaboration Xr] ( |2Q14| ) (for 
fits of T 353 only). 

We initially assume that the proportionality co¬ 
efficients are 30% for the model or the dust map 
itself, and 1 when using (Tr{l^h). These are ad¬ 
justed later as a part of the iterative procedure to 
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obtain a reduced = 1. The is minimized to 
determine the best-fit coefficients, and 

From the results we derive a dust residual map 


D^,,{l,b)=D{l,b) + 

3 

1 ,, ■ N{Ri){l,b)i + fco ■ Wco{l,b) + A°sc 


, 1=1 


( 3 ) 


The map is then filtered to extract the significant 
positive residuals associated with the DNM. The 
residual map is convolved with a Gaussian kernel 
with standard deviation of 1° in order to suppress 
statistical fluctuations while retaining large-scale 
real features. Then we build a histogram of the 
residuals and we fit a Gaussian curve to the core 
of the distribution, within ±2 standard deviations 
of the distribution itself. We set a threshold at 2 
standard deviations of the fitted Gaussian curve, 
and produce the filtered residual map by re¬ 
taining pixels of the original residual map when 
the values in the convolved residual map are above 
the threshold, and setting other pixels to 0. 

Then we proceed iteratively to refine the model. 
At each iteration a = 1, 2,...: 


• The model is defined: 


3 

1=1 

+ ydo-Wco{l,b) + 

+ + ( 4 ) 


• We fit the model to the dust map D 
minimizing the as defined in Eq. 

• The dust residual map is built: 


DZ,{l,b)=Dil,b) + 

■ 3 

ygi,, ■ Ar(Hi)(;, b)i + ^go ■ Wco{l, b) + 


1=1 


( 5 ) 


where and are the best-fit coeffi¬ 

cients; 

• We filter to extract the significant pos¬ 
itive residuals 5^g following the procedure 
described for iteration 0. 


The iterative procedure is stopped when the two 
following conditions are simultaneously satisfied: 


1) the differences between all the best-fit coeffi¬ 
cients in iteration a and iteration a — 1 are com¬ 
patible with 0 within la statistical uncertainties 
from the fit; 2) the difference between minimum 
in iteration a and iteration a — 1 is < 1. 


After the first pass just described, the scal¬ 
ing coefficients for the uncertainties a{l,b) are re¬ 
defined, so that the reduced x^ = 1 the latest 
iteration. Requiring the reduced x^ fo equal 1 
ensures that the statistical uncertainties from the 
fit are truly representative of the differences be¬ 
tween dust maps and models. This, in turn, also 
makes the criteria to stop the iteration of the fit¬ 
ting procedure more meaningful, even if only in a 
data-driven fashion. Using the new scaling coeffi¬ 
cients we perform a second pass, after which the 
final reduced x^ values are found to be 1 within 
1% in all ROIs and variants of the maps. The dis¬ 
tribution of the residuals in the final iteration does 
not show any significant negative values, which 
demonstrates how the procedure is appropriately 
addressing the initial bias due to the missing tem¬ 
plate for the DNM (see Appendix for further 
details). The scaling coefficients for the errors 
on r 353 provided by the |P/ancA; Collaboration XI 
( 2Q14| ) are found to be systematically greater than 
1. The error map represents only statistical uncer¬ 
tainties from a fit where the dust thermal emis¬ 
sion spectral energy distribution was fitted using 
a modified blackbody spectrum. Hence, a scaling 
coefficient greater than 1 may capture the presence 
of systematic errors, or discrepancies between the 
distribution in the sky of dust optical depth and 
our simple model based on the assumption that 
the optical depth scales with gas column densities. 


Because the Planck Collaboration XI (2014) 


recommends the use of R as the best dust column- 
density tracer at high latitudes, we take as our 
baseline DNM templates those obtained from R 
assuming that the error is proportional to the 
model. These maps are shown in Figure We 
note that the DNM component is particularly 
abundant in ROI A, where it follows a known 
interstellar structure dubbed the North Celestial 
Pole (NCP) Loop ( Meyerdierks et al.|1991 ), possi¬ 
bly related to the Ursa Major molecular complex 
(e.g., [Pound fc Goodman 1997). Further results 
from the analysis described in this section are dis¬ 
cussed in Appendix [A} 
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Fig. 6.— DNM templates built from radiance i?, as determined from the iterative fitting procedure for the 
three ROIs A, B, and C. The black lines show the borders of the ROIs. The errors to define the foi* fhe 
fit were assumed to be proportional to the model in Eq. with the proportionality constant defined so that 
the reduced fhe final fit is 1. The maps are shown in zenithal equidistant projection. 
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5. 7-RAY analysis 


5.1. Model Building 


The 7 -ray emission measured using the LAT 
can be modeled as the sum of diffuse compo¬ 
nents and discrete sources. Diffuse emission is 
produced in the interstellar space of the Milky 
Way by energetic interactions of CRs with gas, 
via nuclear collisions leading to 7 -ray production 
mainly through neutral pion decay and by elec¬ 
tron Bremsstrahlung, and with low-energy radia¬ 
tion fields, due to IC scattering by electrons. 


Multiple studies of interstellar 7 -ray emission 
have shown that CR fluxes can be approximated as 
uniform on the scales of interstellar complexes lo¬ 


cated away from potential CR sources (e.g., Abdo 

2QTq1 


et al. 


Planck and Fermi Collaborations 


2014). Under this hypothesis, the 7 -ray intensi¬ 


ties from interactions with gas are proportional to 
the target gas column densities. We can there¬ 
fore model the diffuse intensity as a combination 
of the gas tracer maps described in Section We 
assume the spectral shape to be that measured 
using LAT data for the local interstellar space 
( Casandjian|2Q14' ) in order to perform the fit over 
a broad energy range. This assumption will be 
critically evaluated against the data in the assess¬ 
ment of systematic uncertainties (Sec. 5.4). For 
each tracer map we allow a free scaling coefficient 
in the model, to define the overall CR flux in that 
componenf[^ 

The IC component is modeled based on the 
work by de Palma et al. (2013). They considered 8 
models of IC emission calculated using the GAL- 
PRO fH cod e ( Moskalenko fc Strong|[ 2 QQQ Strong 


et al. 2000 


Porter et al.||2008|) for different val- 
and tuned them 
The work 


ues of some model parameters 
to the LAT data over the whole sky 


by de Palma et al. (2013), originally based on the 
LAT P7 data, was recently updated for a P7REP 
data selection consistent with our dataset. For 
the IC component we allow in the fit a free scal¬ 
ing coefficient, as well as a log-parabolic spectral 
correction to mitigate the effect of uncertainties in 


®HVCs and IVCs span volumes much larger than those of 
interstellar complexes studied so far using LAT data. Were 
there any variations of CR densities within the clouds, the 
scaling coefficient would represent an average CR flux over 
it. 

Tttp://galprop.stanford.edu 


the local electron spectra. Among the models in 


de Palma et al. (2013), we arbitrarily select that 


with the CR source distribution based on the work 
by 


Lorimer 


(2004), a maximum height of the CR 
confinement halo of 10 kpc and a H i spin temper¬ 
ature of 150 K. The other models presented in [d^ 


Palma et al. (2013) will be used later in Section 5.4 


to assess the impact of this choice on the results. 

Finally, there is a diffuse component with al¬ 
most isotropic distribution over the sky inter¬ 
preted as the superposition of extragalactic diffuse 
7 -ray emission and residual backgrounds from CR 
interactions in the LAT misclassified as 7 rays. 
This is modeled using a spectral shape derived 
from the LAT data (see de Palma et al. 2013 for 
the details of the procedure) consistently with our 
IC models, separately for 7 rays converting in the 
front and back sections of the LAT tracker, which 
are characterized by different residual background 
rates relative to their acceptances. Owing to lim¬ 
itations in the assumption that the residual CR 
background can be treated as an isotropic 7 -ray 
term, we allow a free scaling coefficient in the fit 
also for this term. Note that, due to the limited 
size of the ROIs and the quasi-isotropic distribu¬ 
tion of IC 7 -ray emission in the models, both IC 
and isotropic components have to be considered as 
effective background templates in our fits and the 
respective scaling factors are not the object of any 
interpretation in this work. 

We incorporate point sources from the 3FGL 
list {Fermi-LAT collaboration 2015). For each 
ROI we include in the model all the sources within 
10° from the ROI border. The source positions are 
fixed at their catalog values. We assume for each 
source the spectral shape given in the catalog as 
well. For sources outside the ROI but within 10° 
from the ROI border all the spectral parameters 
are fixed at their catalog values. For sources in¬ 
side the ROI we leave all the spectral parameters 
free if the average significance S of the source in 
3FGL is > ^ 1 , just the global normalization free 
if O 2 < S < Oi, and all the parameters fixed oth¬ 
erwise. The thresholds 0i and O 2 are chosen for 
each ROI so that the number of free parameters 
associated with sources is < 35 to prevent con¬ 
vergence problems in the fitting procedure. The 
pairs of thresholds ( 6 >i, 6 > 2 ) chosen are (15cr, 9cr), 
(30cr, 13cr), and (lOcr, 5cr), for regions A, B, and G, 
respectively. We verified that none of the 3FGL 
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sources are positionally coincident with structures 
in the A^(Hi) maps in Figures and|^ Follow¬ 
ing a visual inspection of the residual maps after 
a preliminary analysis we left free in addition to 
those sources also 3FGL J0809.5+5342 in ROI A 
and 3FGL J1159.2+2914 in ROI G. We will show 
later in Figures and that additional two 

years of data in our analysis with respect to 3FGL 
do not introduce any additional sources that sig¬ 
nificantly affect the results of the analysis. 

To summarize, for each ROI the model for the 
7 -ray intensities is described by Equation 




3 

+ ^co • b) qLis{E) + 

+ ^DNM • D{1, b) qLis{E) + 

/ \ <^ic+/3ic In 

+ xic • IC(/, 6, E) f ^ J + 

+ ^iso • Iiso{E) + SRC (6) 


There, gLis(^) stands for the tabulated emissiv- 
ity spectrum of H i in the local interstellar space 
from [Casandjian ( 2Q14[ ). are the Hi col¬ 

umn density maps that we constructed in Sec¬ 
tion |4.1[ so that, for each complex xm is the 
scaling factor with respect to the local emissivity. 
Kinematically-separated complexes along the line 
of sight are numbered as 1 ) low-velocity, i.e., local 
gas, 2 ) and 3) IVCs and HVCs. Wco is the CO 
intensity map described in Section |3.2| (used only 
for ROI A); hence xqo can be interpreted as the 
scaling factor that transforms the local Hi emis¬ 
sivity into 7 -ray emissivity per IFco intensity unit 
in the region studied. Similarly, D is one of the 
templates of the DNM constructed in Section |T^ 
and xdnm is the scaling factor that transforms the 
local Hi emissivity into 7 -ray emissivity per dust 
radiance (or optical depth) unit. IC represents the 
IC emission model: xic is an overall global scaling 
factor, aic and /3ic are the parameters of the log- 
parabolic spectral correction, while Eq is a refer¬ 
ence energy which was set to 700 MeV (the choice 
is arbitrary because the parameter is degenerate 
with the others). Finally, /iso is the tabulated 
isotropic background spectrum (different for front- 
and back-converting events), with xi^o as common 
scaling factor, and SRC represents the model for 
individual sources described above. 


5.2. Analysis Procedure and Baseline Re¬ 
sults 

Using the dataset described in Section [3T| for 
each ROI in Tabled] we build binned count cubes 
with 0?25 spacing in angle and 15 logarithmically 
spaced energy bins from 300 MeV to 10 GeV for 
front- and back-converting events. The model in 
Eq. is fit to the data by using the binned likeli¬ 
hood analysis procedure implemented in the Sci¬ 
ence Tools that takes into account the LAT ex¬ 
posure and PSE. We employ a joint likelihood 
technique to independently treat the front- and 
back-converting events in the fit in order to ex¬ 
ploit the better angular resolution of the former. 
We use the P7REP_V15_CLEAN LAT instrument re¬ 
sponse functions. 

The maximum likelihood values of the coeffi¬ 
cients of the interstellar emission models for the 
baseline analysis are reported in Table In Eig- 
ures[^|^ andj^we show the count maps and resid¬ 
ual maps summed over the entire energy range. 
The residual maps do not show any large-scale 
structures correlated with the templates used to 
model the 7 -ray emission. Some localized residu¬ 
als may hint at missing point sources not included 
in 3EGL, which was developed using the first 48 
months of LAT observations, as opposed to the 73 
months used in our analysis. 

Eor each target (HVG or IVG) in Table we 
calculate the test statistic TS, defined as 

TS = 21„(A), ,7) 


where jC is the maximum likelihood for the model 
that includes the source of interest, while jCq is the 
maximum likelihood of the minimal model where 
the contribution from the source is set to null. TS 
is expected to be distributed as a with a num¬ 
ber of degrees of freedom given by the difference 
of the number of free parameters in the two mod¬ 
els (e.g., Mattox et al.|1996 ). We consider a target 
cloud as detected if TS > 25, which formally corre¬ 
sponds to a 5 cr confidence level that the likelihood 
improvement is not due to statistical fluctuations. 
As detailed in Protassov et al. ( 2002 ), the distribu¬ 
tion of TS may deviate from a if the minimal 
model corresponds to a point on the topological 
border of the phase space of the parameters of the 
model including the source, which applies to our 
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Fig. 7.— Counts (left) and residuals after subtraction of the best-fit baseline model (right) for ROI A. Counts 
and residuals are summed over the entire energy range from 300 MeV to 10 GeV. Residuals are expressed 
in units of approximate standard deviations, calculated as counts minus model-predicted counts divided by 
the square root of counts. The maps are shown in the plate carree projection centered at the ROI center 
used for the y-ray analysis, and are smoothed for display using a Gaussian kernel of cr = 0?5. 
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Fig. 8.— Gounts (left) and residuals after subtraction of the best-fit baseline model (right) for ROI B. See 
the caption of Fig. [^for details. 
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Fig. 9.— Counts (left) and residuals after subtraction of the best-fit baseline model (right) for ROI C. See 
the caption of Fig. [^for details. 


Table 2 

Interstellar best-eit coeeeicients erom the baseline analysis 


Parameter 

ROI A 

ROI B 

ROI C 

XU.1,1 

0.999 ±0.017 

0.94 ±0.07 

1.08 ± 0.04 


0.94 ±0.16 

1.08 ±0.04 

0.68 ± 0.06 

xm,2, 

o o 

d d 
+ 1 

O 

q 

d 

q 

d 

o 

1 + 
o o 

d 

o 


xco"" 

2.3 ± 0.2 

0 .00lg;“ 


a^DNM^ 

0.24 ±0.03 

0.34 ±0.09 

XlC 

1.19 ± 0.05 

1.31 ± 0.05 

1.82 ± 0.09 

otic 

-0.158 ± 0.014 

-0.06 ± 0.02 

0.02 ± 0.05 

^IC 

0.045 ± 0.008 

0.02 ± 0.01 

-0.026 ± 0.014 

^Iso 

0.978 ± 0.017 

1.03 ± 0.02 

0.96 ± 0.04 


Note. —See Eq. I^for a definition of the model parameters. The 
numbering scheme for the xu. 1 ,% coefficients is: 1) low-velocity, i.e., 
local gas, 2) and 3) IVCs and HVCs. 

Note. —Asymmetric error bands are computed using Minos when 
the best-fit value is at the boundary imposed during the likelihood 
optimization process. 

^10^° cm-2 (K km s“^)“^ 
sr 
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situation. However, the TS values for our detected 
targets all greatly exceed the threshold. 

The TS values are reported in Table For 
target HVCs and IVCs that are not detected we 
calculate a 95% confidence level (c.l.) upper limit 
on xhi based on a Bayesian integration of the like¬ 
lihood profile as applied for measurements in the 
2FGL catalog ( Nolan et al.||2Ql^ 3.5). 

In Fig. [^we show the 7 -ray emission compo¬ 
nents associated with detected IVCs, evaluated by 
subtracting from the total 7 -ray counts the com¬ 
ponents attributed to individual sources and fore¬ 
ground/background interstellar emission by the 
likelihood fit. Contours from the V(Hi) maps in¬ 
ferred from the 21 -cm observations in Figures 

and are overlaid. The morphologies of the 7 - 
ray detected signals follow closely the known Hi 
structures in the low-latitude IV Arch and lower 
IV Arch. In the case of the IV Spur a cluster of 
7 rays is associated with the densest part of the 
cloud, but the visual comparison is hampered by 
limited statistics in the 7 -ray dataset. 


5.3. Spectra of the Detected Clouds 


A key assumption of our model is that the emis¬ 
sion of the gaseous components has a spectrum 
described by the emissivity spectrum of local gas 
gnis from Casandjian (2014), modified only by a 
global scaling factor over the whole energy range 
from 300 MeV to 10 CeV. While this assumption 
is expected to hold for the local (low-velocity) gas, 
there may be spectral variations for clouds at great 
distance from the disk. 


For clouds detected in our analysis we can use 
the data to verify how well the ^lis spectrum re¬ 
produces the observed one. To do so we perform 
the 7 -ray analysis independently in 6 broad en¬ 
ergy bins logarithmically spaced from 300 MeV 
to 10 CeV. In these fits the spectral parameters 
are frozen to their determination over the entire 
energy range except for the normalizations. Fig¬ 
ure shows as a reference the scaling factors of 
local Hi emissivity in the three regions studied. 
Figureshows the scaling factors for the emissiv¬ 
ity of all the detected IVCs with respect to local. 





u 





energy (MeV) 


Fig. 11 .— Scaling factor of local Hi emissivity, 
^Hi,i 5 as a function of energy for the three regions 
A, B and C. The shaded bands indicate the values 
obtained from the analysis over the entire energy 
range with their Icr statistical uncertainties. 


The scaling factors are all consistent within sta¬ 
tistical uncertainties with the value determined 
from the fit over the entire energy range. In region 
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Fig. 10.— 7 -ray emission associated with the detected targets: total counts minus the contributions at¬ 
tributed by the fit to all emission components in Eq. but A) the low-latitude IV Arch, B) the lower 
IV Arch, C) the IV Spur. The maps were smoothed using a 2° Gaussian kernel to suppress statistical 
fluctuations and enhance extended features. The overlaid contours of Hi column density have levels A) 
0.6 X 10^^ atoms cm“^ and 1.2 x 10^^ atoms cm“^, B) 0.8 x 10^^ atoms cm“^ and 1.6 x 10^^ atoms cm“^, 
C) 0.9 X 10^^ atoms cm“^ and 1.8 x 10^^ atoms cm“^. 
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Table 3 

Target detections and upper limits 


Region 

Complex 

TS 

95% c.l. upper limit 

A 

Low-Latitude IV Arch 

127 



Complex A 

0 

0.21 

B 

Lower IV Arch 

750 



Upper IV Arch 

0 

0.23 

C 

IV Spur 

154 



Note.— See Eq.|^for the definition of TS. 


Note.— The upper limits are expressed as a fraction of the emis- 
sivity of local H i. 


A the residuals show a negative trend that indi¬ 
cates a spectrum softer than the average in the 
local ISM. Similar trends may be present for both 
the local ISM and IV gas in region B. It is not clear 
if this is real or due to any small imperfections in 
the modeling, but further investigation on this as¬ 
pect is beyond the scope of this work and left for 
future studies. For the scope of our analysis it is 
justified to use the global scaling factors over the 
whole energy range because they are compatible 
within Icr with the weighted average of the values 
for the separate energy bins. 


5.4. Assessment of Systematic Uncertain¬ 
ties 

Uncertainties of inputs to the ^-ray inter¬ 
stellar emission model 

Some inputs to the model in Eq.j^are subject to 
important uncertainties, most notably the DNM 
column densities, IC emission, and the spectrum 
of undetected clouds. We explored the associated 
systematic uncertainties by changing the model¬ 
ing assumptions and evaluating the impact on the 
best-fit parameters for the 7 -ray analysis. 

For each region we modified the baseline DNM 
template by using the alternative templates from 
Section [T^ based on different assumptions for the 
errors in the determination of the DNM template, 
and replacing dust radiance with optical depth. 

We also changed the IC model by using alterna¬ 


tive models as described in de Palma et al. (2013), 


which included assuming maximum heights of the 
CR confinement halo of 4 kpc and 10 kpc, and re¬ 
placing the CR source distribution inferred from 


pulsar observations (Lorimer 2004) with that by 
Case fc BhattacharyaJ^ ( 1998) based on supernova 
remnant observation^ These are two of the most 
important input parameters to GALPROP to de¬ 
termine the morphology and spectrum of IC emis¬ 
sion ( Ackermann et ar]| 2012 a ). 


Finally, we have shown in Section 5.3 that for 
detected targets the assumption that the emission 
spectrum follows the emissivity in the local ISM is 
appropriate, but we cannot do the same for targets 
that are not significantly detected, i.e.. Complex 
A and the Upper IV Arch. We explored possi¬ 
ble variations with respect to the local emissivity 
spectrum based on the set of propagation mod¬ 
els in Ackermann et al. ( 2012 a). We calculated 


the ratios of the emissivity spectra in those mod¬ 
els at all locations in the outer Galaxy (where our 
clouds are located) over the emissivity spectrum 
at the solar circle. We parametrized the extrema 
of these ratios with two functions of energy 


fi = 


E 


1 GeV 


- 0.1 


and 


f2 = 


E 




-[0.08+0.071n(£;/l GeV)] 


( 8 ) 


(9) 


^IGeVy 

Eq. captures the softening due to energy- 
dependent escape that is expected near the bound- 


®We note that the de rivation of the Galactic super nova rem¬ 
nant distribution in |Case Sz Bhattacharya| ( |1998| > is subject 
to uncertain ties and discordant with some later works (e.g., 
|Green|2Q14| >. In our study, though, we use it only as a way 
to probe the uncertainties due to the modeling of CR prop¬ 
agation relying on the previous work by|Ackermann et al.| 
(2012a|);' ^ ' ' 


de Palma et al. 


(20131. 
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energy (MeV) 


Fig. 12.— Scaling factors of H i emissivity in three 
IVCs, xhi ,2 5 as a function of energy in regions A, 
B and C. The shaded bands indicate the values 
obtained from the analysis over the entire energy 
range with their Icr statistical uncertainties. 


aries of the confinement region in the outer Galaxy 
for sources located mostly in the inner Galaxy. In 
addition, Eq. captures a low-energy hardening 
predicted by high-reacceleration models in |Ack 
ermann et al.| (|2012a|) at larger Galactocentric 
mdiET Note that both functions indicate a softer 
spectrum. In order to be conservative we explore 
the extreme assumptions that the emissivity spec¬ 
trum may actually be harder at the locations of 
the undetected targets by inverting the signs of 
the exponents 


- (nsv) ■ <“> 

and 

+ [0.08+0.07 ln(£;/l GeV)] 

• ( 11 ) 

We use these functions to modify the emissivity 

spectra of the undetected clouds, by multiplying 
the local emissivity spectrum g^Lis in Eq- by the 
correction function and scaling it so that the emis¬ 
sivity integrated in the 300 MeV to 10 GeV range 
is the same as for gLis- In this way, we can still in¬ 
terpret the upper limits on xhi for the undetected 
components as scaling factors with respect to the 
local emissivity, and, at the same time, explore the 
impacts on the other fit parameters. 

We repeat the baseline analysis described in 
Section [5^ for all the combinations of the model 
variations described above, and, for each param¬ 
eter, we take the extreme variations with respect 
to the baseline fit as our estimate of the system¬ 
atic uncertainty related to the model. In the cases 
of upper limits for undetected targets we use the 
worst, i.e., largest, upper limits in the discussion 
that follows. 

We summarize the scaling factors for H i emis¬ 
sivity in the different regions and complexes in Ta¬ 
ble including the spread due to inputs to the 
7 -ray interstellar emission model as determined 
through the methodology described in this sub¬ 
section. In general, the emissivity spread for the 



^ Our choice to include this effect in the evaluation of 
systematic uncertainties is conservative, because the low- 
energy hardening is a feature that appears only in mod¬ 
els with high reacceleration, while models with convection 
or pure diffusion models are also viable (| Ackermann et al.| 
|2012a| >. 
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detected complexes is driven by the choice of the 
DNM template, and, to a lesser extent, the choice 
of IC model for local gas (which is less structured, 
hence more degenerate with IC emission). For 
non-detected complexes the systematic uncertain¬ 
ties are dominated by the choice of the emissivity 
spectrum. The results in ROI B show a stronger 
dependence on the assumed IC model, which af¬ 
fects also the determination of the lower IV Arch 
emissivity and the upper limit on the emissivity of 
the upper IV Arch at 2 statistical a level. 


5.4-2. Jackknife Tests 


So far we have characterized the statistical un¬ 
certainties in the fit parameters in Section 5.2 and 
systematic uncertainties related to some inputs to 

However, other poten¬ 


the model in Section 5.4.1 


tial sources of systematic uncertainties are related 
to the assumption that the emissivity scaling co¬ 
efficients for each component are constant across 
a given ROI, and to the presence of regions with 
large deviations or mismodeled individual sources 
driving the fits. 

To characterize the magnitude of these uncer¬ 
tainties we perform jackknife tests where we repeat 
the analysis described in Section [5^ 150 times for 
each ROI, each time masking a circular region with 
random center within the ROI and covering 20% 
of the ROI area. The mask size was chosen to be 
much larger than the 68% event containment ra¬ 
dius of the LAT for front-converting photons at 
the lowest end of the energy range at 300 MeV 
1?5) and large enough to probe possible mild 
variations of the CR densities across the complexes 
studied. 


Figures p!3lp^ and show the distributions of 
the HI emissivity scaling coefficients for the three 
ROIs. The distributions all peak within la from 
the best-fit values obtained from the baseline anal¬ 
ysis and have standard deviations smaller than or 
comparable to the la statistical uncertainties from 
the 7 -ray fits. This shows that from a statistical 
point of view the obtained parameter values are 
representative of the entire ROIs. 

In ROI A, Figure the distribution of the 
scaling coefficients for Complex A, xhi, 3 , narrowly 
peaks at zero, showing that this parameter is al¬ 
ways at the boundary of the allowed values. In 
the rest of the work, we will make use only of the 


upper limit on the Complex A emissivity that, on 
the contrary, is a robustly determined quantity. 
Similar considerations hold for the upper IV Arch 
in ROI B. In all ROIs, most notably in ROI A, 
the distribution of the scaling coefficients xhi ,2 for 
the IVCs show tails at low or high values, which 
suggests that emissivities lower or higher than av¬ 
erage are possible within the gas encompassed by 
the respective templates. 

To take into account these uncertainties in our 
estimates of the emissivities we compute, based 
on the distributions of the fit parameters in Fig- 
ures pAj andpls) the la confidence intervals for 
detected complexes and 95% c.l. upper limits for 
complexes that are not significantly detected, as 
reported in Table 0 (details of the derivation are 
given in Appendix]^ . For the detected complexes 
the confidence intervals are of the same magnitude 
as the statistical errors. Hence we will combine 
them with the other uncertainties for the discus¬ 
sion that follows. In the case of the complexes that 
are not significantly detected the sources of uncer¬ 
tainty probed through the jackknife tests (such as 
possible variations of the CR densities across the 
complexes studied) have an impact on the emis¬ 
sivity upper limits much smaller than those due 
to uncertainties in the analysis model inputs and 
those encoded by the statistical error from likeli¬ 
hood fit, so they will be neglected. 


6. RESULTS AND DISCUSSION 


6.1. The 7 -Ray Emissivities of HVCs and 
IVCs 


The Hi emissivity scaling coefficients in Table 
can be converted into emissivity values. To do 
so we multiply them by the integrated emissivity 
betwe en 300 MeV and 10 GeV from [Casandjian 
(2014), i.e., 7.1 x 10“^^ s“^ sr“^ H“^. Note that 
this rests on the measurement described in Sec¬ 


tion 5.3 that the shape of the emissivity spectrum 
is the same as in the local ISM for detected com¬ 
plexes, and is valid within the extreme variations 
parametrized in Eq. for those that are 

not significantly detected. 


The conversion into emissivities requires tak¬ 
ing into account one additional source of sys¬ 
tematic uncertainty, i.e., the characterization of 
the LAT instrument response, and in particu¬ 
lar its effective area. Based on the estimate of 
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^HI,3 


Fig. 13.— Distribution of the xhi scaling coefficients in ROI A from the jackknife tests. In each panel the 
solid vertical line marks the best-fit value from the baseline analysis, and the dashed lines delimit its ±lcr 
statistical uncertainty interval. 
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^HI,3 


Fig. 14.— Distribution of the xm scaling coefficients in ROI B from the jackknife tests. See caption of 
Fig. for details. 
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Table 4 

Hi Emissivity scaling factors 


Region 

Complex 

Scaling Factor 

Stat^ 

Sys (Model)^ 

Sys (Jackknife)*^ 

A 

Local 

Low-Latitude IV Arch 
Complex A 

0.99 

0.94 

< 0.21 

0.02 

0.15 

+ 0.04 
- 0.14 
+ 0.26 
- 0.01 
< 0.22 

0.08 

0.08 

< 0.001 

B 

Local 

0.94 

0.07 

+0 . 13 

- 0.22 

0.08 


Lower IV Arch 

1.08 

0.04 

+ 0.10 

- 0.08 

0.09 


Upper IV Arch 

< 0.23 


< 0.45 

< 0.01 

C 

Local 

1.08 

0.04 

+ 0.05 

- 0.09 

0.09 


IV Spur 

0.68 

0.06 

+ 0.09 
- 0.03 

0.06 


Note. —Upper limits are provided at 95% confidence level. 

^Icr statistical uncertainty from the likelihood analysis. 

^ System atic spread from varying some inputs to the 7 -ray interstellar emission model (see 
Sec. |'5.4.1| ). 

‘^Icr u ncertai nty or 95% c.l. upper limit from the distributions obtained with the jackknife test 
(see Sec. |5.4.'^ . 




Fig. 15.— Distribution of the xhi scaling coefficients in ROI C from the jackknife tests. See caption of 
Fig. for details. 
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the systematic uncertainties of the LAT effective 
area in Ackermann et al. ( 2012 c), the system¬ 


atic error on the local emissivity by |Casandjian| 
(2014) integrated between 300 MeV and 10 GeV 


-27 ^-1 jj-1 


is 0.6 X 10 

Note also that in many HVCs and IVCs the ra¬ 
tio of the ionized gas mass to the neutral atomic 
gas mass is much larger than in the local ISM 


and close to unity (e.g., Wakker et al. 2008). If 
the column densities of ionized gas are correlated 
with those of neutral gas the scaling factors in 
the 7 -ray fits will be overestimated to encom¬ 
pass the emission from ionized gas. In this re¬ 
spect our estimates of scaling factors and emis- 
sivities should be regarded as an upper limit to 
the real values. The same consideration applies to 
the potential presence of CO-dark molecular gas 
in HVCs/IVCs. Aside from the emissivity scaling 
coefficient slightly larger than 1 for the lower IV 
Arch there is no clear indication of such effect from 
our analysis. 

The final results for the emissivities of the tar¬ 
gets are summarized in Table For the purpose 
of interpretation we will rely on those findings and 
on the emissivity scaling factors in Table which 
to a good approximation are not subject to un¬ 
certainties due to the LAT instrument effective 
area. 

The 7 -ray emissivities are a proxy for the 
CR densities in the complexes studied. Ow¬ 
ing to the cross sections for 7 -ray production 


in nucleon-nucleon inelastic collisions (e.g., Ka- 
mae et al. 2006), the 7 -ray energies studied be¬ 


tween 300 MeV and 10 GeV constrain the CR 
densities from energies of ^ 3 GeV/nucleon up 
to ^ 200 GeV/nucleon. Note that, while we do 
not expect significant changes in the He/H ratio 
because both elements are primordial, the abun¬ 
dances of metals, i.e., elements heavier than He, 
are known to vary in HVCs and IVCs, and depend¬ 
ing on the metallicity of target gas the emissivity 
per H atom can vary up to < 5% of the total 


(e.g., Mori 2009). Furthermore, the total emis¬ 
sivity depends also on the CR composition and 
spectra of the different elements at a level up to 
20% to 30% of the total emissivity per H atom 


(Kachelriess et al.||2014). Part of the emissivity 




variations (Table or, equivalently, differences in 
the emissivity scaling factors (Table therefore 
can be attributed to variations in the CR elemen¬ 


tal composition and of their spectra, as well as to 


the uncertain metallicity of the targets (Wakker 


2001 ). 


Our measurement of the local gas emissivity in 
each ROI agrees well with unity, indicating the 
robustness of the analysis procedure and assump¬ 
tions. The emissivity of the IV Arch indicates a 
^ 50% decrease of the CR densities with respect 
to the local value within 2 kpc from the Galac¬ 
tic plane, while the HVC Complex A provides the 
strongest limit on the CR densities at a few kpc 
above the disk. The implications of these results 
are discussed in the following sections. 


6.2. Testing the Origin of Cosmic Rays in 
the Galactic Disk 

If CRs originate in the disk of the Milky Way 
and diffusively propagate in the surrounding halo 
we expect the CR densities, and hence the 7 -ray 
emissivities, to decrease with distance 2 ; from the 
disk itself^ We test this hypothesis against our 
data in Table by computing the Kendall r rank 
correlation coefficient. Using the emissivity scal¬ 
ing coefficients in Table and adopting for z the 
mean within the bracket for targets in Table 
and 2 ; = 0 kpc for the local gas in each ROI, we 
obtain r = —0.54. The negative value indicates, 
indeed, a negative correlation of emissivity with 

ki- 

We calculate the significance of this trend by 
comparing the r of the actual data with a dis¬ 
tribution of the r coefficients obtained from the 
null hypothesis of no correlation between emissiv¬ 
ity and 2 :. We generate null hypothesis datasets 
starting from the real dataset using the following 
procedure: 

• we draw a set of emissivity scaling coeffi¬ 
cients with normal distribution centered at 
the measured values and with a equal to the 
sum in quadrature of statistical uncertain¬ 
ties and systematic uncertainties from the 
jackknife tests in Tablefor significantly de¬ 
tected targets and local foreground gas in 
each ROI; 


Owing to the distances to other massive star-forming galax¬ 
ies, any contributions from CRs escaping from even the 
closest systems would be too small to affect this conclu¬ 
sion and also smaller that the uncertainty of our emissivity 
measurements. 
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Table 5 

Emissivities of HVCs and IVCs 


Region 

Complex 

z 

(kpc) 

7 -ray emissivity 

Stat^ 

Sys (LAT)‘> 
sr“ 

Sys (model)*^ 
H-i) 

Sys (Jackknife)*^ 

A 

Low-Latitude IV Arch 
Complex A 

0 . 6 - 1.2 

2 . 6 - 6.8 

6.7 

1.1 

0.6 
< 1.7 

+1.8 

-0.1 

0.6 

B 

Lower IV Arch 
Upper IV Arch 

0.4 1.7 

0.7-1.7 

7.7 

0.3 

0.7 

< 3.5 

+0.7 

-0.6 

0.6 

C 

IV Spur 

0.3-2.1 

4.8 

0.4 

0.4 

-0.2 

0.4 


Note. — 2 : brackets are replicated from Table ^ for the reader’s convenience. Data sources are listed there. 

Note. — 7 -ray emissivities are integrated between 300 MeV and 10 GeV. 

Note. —We provide 95% confidence level upper limits for the worst-case scenario from the model variations considered and 
taking into account uncertainties due to the LAT effective area. 

^Icj statistical uncertainty from the likelihood analysis. 

^Error from uncertainties in the LAT effective area. 

‘^Systematic spread from varying some inputs to the 7 -ray interstellar emission model (see Sec. |5.4| ). 

‘^Icr uncertainty or 95% c.l. upper limit from the distributions obtained with the jackknife test (see Sec.jsT^. 


• we add upper limits at the values measured 
in Table for targets which are not signifi¬ 
cantly detected; 

• we randomly shift both scaling factors and 
upper limits with uniform probability dis¬ 
tribution within the systematic uncertainty 
brackets from the model input variations in 
Table (assuming the bracket to be sym¬ 
metric for upper limits); 

• we draw a set of z values with uniform proba¬ 
bility distribution within the brackets in Ta¬ 
ble and add three elements with z = 0 
to represent the local foreground gas in each 
ROI; 

• we randomly shuffle in an independent fash¬ 
ion the emissivity scaling factors and z sets 
so obtained to be consistent with the null 
hypothesis of no correlation between the two 
quantities. 

The distribution of 30000 null hypothesis datasets 
shows that there is a 2.5% chance probability to 
obtain a r coefficient < —0.54. 

Therefore, our measurements provide evidence 
at 97.5% confidence level that the 7 -ray emissiv¬ 
ity per H atom decreases with distance from the 


disk of the Galaxy. This corroborates the notion 
that CRs at the relevant energies from GeV to 
TeV originate in the Galactic disk by directly trac¬ 
ing the distribution of GRs in the halo for the 
first time. Our result goes far beyond the test of 
the Galactic origin of GRs proposed by |Ginzburg| 
(1973), i.e., measurement of the 7 -ray emissivity of 


Magellanic Glouds, that was performed by lSreekiT 


mar et al. (1993). 


We note that we also expect a decrease in gas 
metallicity, hence emissivity per H atom, with dis¬ 
tance from the disk if the most distant targets like 
Gomplex A are primeval gas falling into the Milky 
Way ( Wakker|2QQl ). However, the expected effect 
of at most ^ 5% on the 7 -ray emissivities (Mori 


2009) is not sufficient to explain the magnitude of 


the variations observed between the scaling coef¬ 
ficients in Table Ghanges in the spectra of the 
different GR elements could also cause emissivity 
variations. However, the magnitude of this effect 
is estimated to be at most 20% to 30% of the total 
emissivity per H atom (taking into account He in 
both GRs and the ISM), lower than the decrease 
of emissivity of > 50% for the Upper IV Arch and 
of > 75% for Gomplex A that we measured. 

Additionally, in the region of the outer Galaxy 
where our targets are located the distance to the 
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center of the Milky Way increases with increas¬ 
ing distance from the disk as well. While some 
propagation models advocate a sizable decrease 
of CR densities with increasing Galactocentric ra¬ 


dius, observations (Abdo et ah 2010 Ackermann 
et al.|[2MT ) indicate that the gradient of CR den¬ 


sities is marginal up to 15 kpc from the Galactic 
center. We conclude that the observed emissiv- 
ity decrease as a function of ^ is more likely due 
to diffusion in the halo. A quantitative compar¬ 
ison with propagation models follows in the next 
section. 


6.3. Comparison with CR Propagation 
Models 

We compare our results with predictions from 
GALPROR We will rely on the set of GALPROP 


models in Ackermann et al. ( 2012 a) that have been 


compared extensively to direct CR measurements 
and LAT data, and which include a number of 
variations of the most relevant and uncertain pa¬ 
rameters of the calculations. In I Ackermann et al.l 


( 2012 a) the CR transport equations were solved 


for the case of cylindrical symmetry, based on the 
assumption that CRs are injected in the disk of 
the Milky Way and propagate in a cylindrical vol¬ 
ume with boundaries at maximum radius from the 
center of the Galaxy Rmax and maximum height 
from the disk ^^max where the CR particle number 
densities are required to go to zero. 

Figure shows the emissivity scaling factors 
fr om Table against predictions from the models 
in Ackermann et al. ( 2012 a). The model predic¬ 
tions are calculated for distances along the line of 
sight approximately corresponding to the peak gas 
column densities for each ROf^ We normalized 
the model emissivities to the value in the disk at 
the solar circle. The model emissivities at the solar 
circle integrated from 300 MeV to 10 GeV in 7 - 
ray energy were found to be within ^ 7% from the 
measurement in Casandjian (2014). Note that the 
calculation of the emissivities in I Ackermann et al.l 


( 2012 a) only accounts for interactions between CR 
and interstellar gas nuclei with atomic number 


^^This accounts for the different ranges of Galactocentric ra¬ 
dius in each line of sight. The effect is clearly visible for 
ROI A in Figurewhere the profiles are concave, and the 
curves for Zmay^ = 10 kpc bifurcate into two families for 
Rmax = 20 kpc and Rmax = 30 kpc due to particle escape 
at the radial boundary. 


< 2. Rescaling to the locally measured emissivity 
takes into account the contribution from metals 
in the targets and heavier CR species. Variations 
of the target gas metallicities are not taken into 
account in these models, and could produce a fur¬ 
ther decrease of emissivity up to at most 5% with 
increasing 2 ; ( Mori|[2QQ9 ). 


The model parameter with the largest impact 
on the vertical gradient of CR densities, or, equiva¬ 
lently, 7 -ray emissivities, is the maximum height of 
the confinement halo ^^max- Models in | Ackermann] 
assume values of Zmax between 4 kpc 
While there is a broad agreement 
between models and measurements, the rapid de¬ 
crease of the emissivities within 2 kpc from the 
disk inferred from the upper limit for the upper 
IV Arch seems to favor smaller values of ^^max- 
The lower and upper IV Arch have similar altitude 
brackets, 0.4 — 1.7 kpc and 0.7 — 1.7 kpc, respec¬ 
tively, but their emissivities are consistent with 
local for the first and < 50% of local for the lat¬ 
ter. This is not necessarily contradictory because 
the distance brackets are pairs of upper and lower 
limits on the distances; hence the two clouds may 
be separated by a physical distance up to ^ 1 kpc. 


et al. ( 2012 a 


and 10 kpc. 


To obtain a crude estimate of from our 
results we fit to the emissivity scaling factors of 
IVCs and HVCs from Table ID a function 


e(z) = 1 




'^^max 


/3 


( 12 ) 


where 5 is the emissivity expressed as a fraction 
of the emissivity at ^ = 0, 2)niax is the maxi¬ 
mum height of the confinement halo and (3 is an 
index which modifies the steepness of the verti¬ 
cal gradient. Eq. with /3 = 1 is an accurate 
^(-dependent part of the solution of the diffusion 
equation in the plane-parallel geometry (infinitely 
thin Galactic plane) with a uniform source dis¬ 
tribution when only ionization losses are assumed 
(e.g., Jones et ar]2QQl|). T he GALPROP models 
in [Ackermann etah (2012a) have vertical emissiv¬ 
ity gradients at given Galactocentric radius char¬ 
acterized by 1 < /d < 1.5. Values of P greater than 
1 result from energy dependent escape at the ra¬ 
dial boundary of the propagation model and, to a 
lesser extent, from energy losses of nuclei. In any 
case, we stress that the functional form in Eq. 
holds only in diffusion-dominated scenarios where 
it is required that the CR densities go to zero at 
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(A) 



(B) 



Fig. 16.— For ROIs A, B, and C we compare the emissivity scaling factors in Table (gray rectangles) 
with predictions from the models in Ackermann et al. (2012a) (curves). The horizontal widths of the 
rectangles indicate the ^ brackets of target IVCs and HVCs, i.e., the range between lower and upper limits 
on their altitudes (Table . The dark gray rectangles have vertical size corresponding to the statistical 
uncertainties, while for the light gray rectangles the vertical size encompasses the sum in quadrature of 
statistical uncertainties and systematic uncertainties from the jackknife tests, augmented by the systematic 
spread from variations to the analysis model inputs from Table The emissivity of local gas is assigned 
to the range from z = 0 kpc to ^ = 0.3 kpc (disk). The model curves from Ackermann et al. (2012a) were 
calculated for the line of sight indicated in the legend of each panel, approximately corresponding to the 
column density peaks of the target complexes. The curves are color-coded based on the maximum heights 
^max of the CR confinement halo in the models. 
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2 ^ ^max' 

We perform the fit by using a maximum likeli¬ 
hood method where we neglect the different Galac- 
tocentric distances of the targets. We assume for 
the emissivity scaling factors a Gaussian probabil¬ 
ity distribution with a given by the quadratic sum 
of statistical uncertainties, jackknife uncertainties, 
and the largest brackets from model input varia¬ 
tions from Table Upper limits including sys¬ 
tematic uncertainties from Table [4] are assumed 
to be the 95% containment value of a Gaussian 
probability distribution. Finally, we assume that 
the probability distribution of target distances is 
uniform within the brackets in Table Owing 
to the paucity of the data points, rather than fit¬ 
ting the index [3 to the data we assume the ex¬ 
treme values P = I and /3 = 1.5 as found from 
the fits to GALPROP predictions. The best fit to 
the data is obtained for = (2.2 ±1.9) kpc, 
and Zmax = (1.8 ±2.1) kpc for the two values of /3. 
Note that the profile would be steeper, hence ^^max 
smaller, if the emissivities were biased by the pres¬ 
ence of unaccounted ionized gas in HVGs/IVGs. 

This method provides an estimate of ^^max < 
6 kpc (2cr upper containment value from the fits 
above). Most of the constraining power on the 
^max parameter, however, at present comes from 
the upper limit on the emissivity of the upper IV 
Arch. The latter hints at a possible tension with 
indications from the modeling of other observables 
related to GRs. Differences with respect to studies 
of synchrotron emission ( Orlando fc Stron^|2Q13 ) 
can be understood if the densities of GR electrons 
beyond the boundaries of the confinement halo 
are sufficient to produce the emission observed in 
the domain from radio to microwaves. However, 
the interpretation of the fiat radial profile of y-ray 
emission produced by interactions of GR nuclei in 
the outer disk of the Milky Way in terms of a 
large ^^max of tho order of ^ 10 kpc in the context 
of GALPROP models similar to those considered 


here (Abdo et al. 2010 Ackermann et al. 2011) 


seems to be disfavored 

The observed emissivities of the HVG and IVGs 
were compared to predictions from the limited 
set of diffusive reacceleration models for a sim¬ 
plified axisymmetric model of the Galaxy in | Ack- 
ermann et al. (2012a). Other models proposed in 


the literature include GR-driven Galactic winds 


and anisotropic diffusion (e.g., Breitschwerdt et al. 


2002), and a non-uniform diffusion coefficient that 


increases with the Galactocentric radius and the 
distance from the Galactic plane (e.g., Shibata 
et al.||2007 ), as well as more sophisticated ways of 


modeling the Milky Way structure (e.g., Werner 
et al. poTsj ), but further comparison to theoretical 


models is beyond the scope of this article. 

Finally, we note that the GR propagation equa¬ 
tions are often solved assuming that the GR par¬ 
ticle number densities go to zero at the bound¬ 
aries of the propagation volume. This assumption 
has little effect on the predicted GR fluxes for the 
Galaxy as a whole and at a considerable distance 
from the boundaries. However, close to the halo 
boundaries one may expect significant deviations 
from the model predictions and tracing the GR 
distribution in the halo can provide unique model- 
independent information about its structure and 
extent. Studies of the GR outflow from the Milky 
Way into the intergalactic space that are still in 
their infancy (e.g., Everett et al.||2Q12 ) could also 
benefit from the method illustrated in this work. 


7. CONCLUSIONS 

We have searched for high-energy y-ray emis¬ 
sion from a sample of HVGs and IVGs in the halo 
of the Milky Way at different distances from the 
Galactic disk using 73 months of data from the 
LAT in the energy range between 300 MeV and 
10 GeV. 

We have achieved the first detections of IVGs 
in y rays, the low-latitude IV Arch, the lower IV 
Arch, and the IV Spur, and set constraints on the 
emission from the remaining targets, the upper IV 
Arch and the HVG Gomplex A (Table [^. The 
spectra of the detected complexes were all consis¬ 
tent within statistical uncertainties with that of 
gas in the local interstellar space. 

We And evidence at 97.5% confidence level that 
the y-ray emissivity per H atom of the clouds de¬ 
creases as a function of distance from the disk. 
This corroborates the notion that GRs are accel¬ 
erated in the Galactic disk and then propagate in 
a surrounding halo, for the first time in a direct 
way. 

The y-ray emissivity per H atom of the upper 
IV Arch hints at a 50% decline of the GR densities 
within 2 kpc from the disk. The upper limit on the 
emissivity of Gomplex A at 22% of the local value 
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gives the most stringent constraint to date on the 
fluxes of CRs at few kpc from the Milky Way disk. 

Our results can be compared to CR propagation 
models and inform further development to more 
realistically take into account the escape of CRs 
from the halo of the Milky Way and the outflow 
of particles merging into the local intergalactic 
medium. At the same time, the observational con¬ 
straints can be improved if distance brackets for 
more HVCs and IVCs in the mass-distance range 
detectable by the LAT are measured or existing 
brackets are made more constraining, e.g., by the 
ongoing Gaia survey (e.g., de Bruijn^|2Q12 ). 
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A. ADDITIONAL RESULTS ON THE INTERSTELLAR MEDIUM 


The analysis described in this paper is aimed at determining the CR content of IVCs and HVCs. Never¬ 
theless, we obtained results relevant to the properties of interstellar gas and dust in such clouds, as well as 
in the local complexes seen in their foregrounds. This appendix summarizes these results, including some 
methodological aspects. 

A.l. Example of Results from the Iterative Dust Fitting Procedure 

In this section we take as an example ROI A to illustrate some details of the iterative fitting procedure 
applied to the dust maps in Section [4^ We follow the fit of the radiance map, R, for the assumption that 
the errors are proportional to the model M from Eq. |^in calculating the The results from all iterations 
are shown in Figures pT| pisl and[^ 

F ignres pT| and \T8\ show how the scaling coefficients of the Hi and CO maps change significantly during the 
iterative procedure and how they stabilize at a final value at the end. This can be interpreted as removing the 
bias due to the missing DNM component in the initial iteration, which is alleviated as the missing component 
is recovered based on the data themselves. 

Figure shows the difference between the filtered dust residual maps from the last iteration and the 
initial one. It is evident that borders exist between the Hi-dominated and DNM-dominated, and between 
the DNM-dominated and CO-dominated regions, where the missing component is not recovered from the 
residuals in the initial iteration. This can be understood because the Hi and CO components in the initial 
fit compensating for the missing DNM component. 

Finally, we note that while in the initial iterations of the fits there are often both positive and negative 
significant residuals, only positive ones are found at the end of the iterative procedure. This is consistent 
with our working hypothesis that the missing DNM component traced by the dust maps is recovered by the 
template determined from the iterative procedure. 


A.2. Fitting Coefficients for the Dust Maps 


The final fit coefficients, reduced values, and error scaling factors for all the fits performed are given 
in Tables miziiiiioi andpr] In ROI B there are no significant positive residuals when fitting T 353 , so the 
iterative procedure was not performed. The number of iterations varies depending on the region and the 
number of free parameters considered, ranging from at least five in ROI B up to at most 30 in ROI A. 


The scaling factors for the errors on r 353 given by the Planck Collaboration XI (2014) are ^ 5. Possible 


Table 6 

Results from R fit in ROI A 



a oc R 

a oc M 

a scaling factor 

15% 

12 % 

Number of iterations 

12 

11 

xVn.d.i. 

1.006 

1.000 

yni,P 

2.160 ±0.008 

2.025 ± 0.007 

yni,2'' 

3.19 ±0.04 

2.41 ±0.03 

yni,3^ 

0.00 ± 0.01 

0.000 ± 0.006 

yco^ 

3.84 ± 0.18 

4.48 ± 0.10 

yres 

1.001 ± 0.002 

1.000 ± 0.007 

F? c 

-fi-iso 

0.731 ± 0.019 

0.811 ± 0.016 


^10-32 sr-i 

W cm“^ sr“^ (K km 
W m-2 sr-^ 
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iteration 


Fig. 17.— Best-fit parameters of the model described in Eq. for ROI A in each iteration of the dust 
analysis. We show as an example the case for which the radiance map R is considered and the errors in the 
definition are assumed to be proportional to M. 
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iteration 


Fig. 18.— Variation of the best-fit parameters of the model described in Eq. |^for ROI A in each iteration 
of the dust analysis with respect to the previous one (expressed as the number of la statistical uncertainties 
from the fits). We show as an example the case for which the radiance map R is considered and the errors 
in the definition are assumed to be proportional to M. 
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Fig. 19.— Normalization of the residual dust component, ^res, (left), and value (right) for ROI A in each 
iteration of the dust analysis. We show as an example the case for which the radiance map R is considered 
and the errors in the definition are assumed to be proportional to M. 
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Fig. 20.— Difference of the filtered dust residual map from the last iteration minus the map from iteration 
zero, divided by the former. We show as an example in ROI A the case for which the radiance map R is 
considered and the errors in the definition of the y^ are assumed to be proportional to M as defined in Eq. 
The map is shown in the plate carree projection used in the dust and y-ray fits for ROI A, and the magenta 
line shows the border of the region considered in the fit. 
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Table 7 

Results from T353 fit in ROI A 


cr OC CTt- CT OC r353 cr OC M 


cr scaling factor 
Number of iterations 

xVn-d-f- 

yco^ 

yres 

T- ^ 

' ISO 


4.9 

11 

0.993 

0.738 ± 0.002 
0.879 ±0.013 
0.01 ± 0.03 
2.30 ±0.08 
0.995 ± 0.014 
0.00 ± 0.01 


23% 

9 

1.003 

0.728 ± 0.002 
0.859 ± 0.016 
0.00 ±0.01 
2.31 ± 0.09 
0.997 ± 0.013 
0.000 ± 0.012 


24% 

5 

1.000 

0.806 ± 0.003 
0.899 ±0.021 
0.000 ±0.006 
3.2 ± 0.2 
1.001 ± 0.001 
0.000 ± 0.005 


^10“2® cm^ 

(K km s“^)“^ 
n0“® mag 


Table 8 

Results from R fit in ROI B 



a OC R 

a OC M 

cr scaling factor 

14% 

17% 

Number of iterations 

6 

5 

xVn.d.f. 

0.996 

1.008 

yni,!"" 

2.187 ± 0.010 

2.246 ± 0.013 

yni,2^ 

1.982 ± 0.011 

2.168 ± 0.014 

ym,3'^ 

0.46 ± 0.03 

1.07 ± 0.04 

yres 

0.986 ±0.015 

0.98 ±0.06 

F? ^ 

-TLiso 

0.935 ± 0.012 

1.039 ± 0.016 


^10-32 sr-^ R-^ 
bio-io W m-2 sr-^ 


Table 9 

Results from T353 fit in ROI B 


cr OC (Tt cr OC r 353 a OC M 


cr scaling factor 
Number of iterations 
xVn-d.f. 

yiii,2^ 

yni,3'" 

yres 

r- 

' ISO 


7.0 

5 

0.994 

0.729 ± 0.005 
0.390 ±0.005 
0.146 ±0.012 
1.03 ± 0.02 
0.50 ±0.06 


36% 

6 

1.003 

0.680 ± 0.003 
0.351 ± 0.003 
0.152 ± 0.012 
0.98 ± 0.02 
0.00 ±0.01 


32% 

8 

0.994 

0.708 ± 0.006 
0.401 ± 0.007 
0.221 ± 0.018 
1.07 ± 0.06 
2.60 ±0.08 


^10“2® cm^ R-i 
’^lO-^ mag 
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Table 10 

Results from R fit in ROI C 



a oc R 

a oc M 

a scaling factor 

11 % 

13% 

Number of iterations 

6 

7 

X^/nM.f. 

1.058 

1.027 

yni,i^ 

2.045 ± 0.011 

2.259 ± 0.013 


2.124 ± 0.019 

2.32 ± 0.02 

yres 

1.000 ±0.014 

1.00 ±0.04 

P b 

-fl'lSO 

1.38 ±0.02 

1.25 ±0.02 


^10-32 sr-1 H-i 

biQ-io w m-2 sr-^ 


Table 11 

Results from T353 fit in ROI C 



cr OC CTt- 

o- oc r353 

a oc M 

a scaling factor 

5.4 

23% 

26% 

Number of iterations 

8 

7 

6 

xVn-d-f- 

1.018 

1.018 

0.994 

yni,i'' 

0.745 ± 0.003 

0.682 ± 0.003 

0.865 ± 0.004 

yiii,2'" 

0.320 ±0.006 

0.300 ±0.005 

0.377 ± 0.007 

yres 

1.001 ± 0.004 

0.992 ± 0.011 

1.033 ±0.039 

Tiso^ 

0.000 ± 0.005 

0.000 ±0.003 

0.000 ±0.006 


^10-26 cm^ 
mag 
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reasons for this are discussed in Section 4^ The scaling factors obtained when assuming a oc R are 10% 
to 15%; when assuming a oc rsss they are ^ 25%. It can be expected that the dispersion on R is smaller 
because this is an integrated quantity rather than referring to a specific frequency. 

We note that the ym,i coefficients can have a physical interpretation in terms of average specific power, 
i.e., power radiated by dust per H atom, per unit solid angle, R/N(R), and average specific opacity, i.e., 
optical depth of the dust per H atom, r 353 /A/'(H). Non-zero coefficients for IVCs are consistent with previous 
measurements of dust thermal emission from this class of objects (e.g., Woerden et al. 2QQ5[ Chapter 
9, and references therein), and, for ROIs A and B, with the measurements of metallicities close to solar 
(Table [^. On the other hand, the ym ,3 coefficient for Complex A being compatible with 0 in Tables and 
is consistent with the low metallicity measured in this cloud ( Wakker||2QQl ). 


A.3. Differences between Alternative Determinations of the DNM Column Density Tem¬ 
plates 


In Figure we show the relative difference between different determinations of the DNM column density 
templates for ROI A, namely the difference between the baseline determination of the DNM column densities 
based on R assuming a (x M (see Eq. I^for a definition of the quantities), and from two alternative determi¬ 
nations based on R assuming a x R (left), and based on r 353 assuming a x M (right). For the second case 
we calculated the average R/r^^s ratio in non-empty pixels of the two maps and used it to convert r 353 in 
equivalent R. 


The figure illustrates how the different assumptions influence the final DNM template, for example how 
some regions have DNM only for one of the cases and not the other, and how the contrast within the map can 
change. For these reasons we considered the different determinations of the DNM template in the assessment 


of the systematic uncertainties for the y-ray analysis in Section 5.4 
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Fig. 21.— Differences between alternative determinations of the DNM column density templates for ROI A. 
Left: difference between the DNM template derived from R assuming a x M minus the DNM template 
derived from R assuming a oc R, divided by the former. Right: difference between the DNM template 
derived from R assuming a x M minus the DNM template derived from r 353 assuming a x M (converted 
into R by multiplying by 1.88 x 10“^ W m“^ sr“^ mag“^, the average of the ratios between non-zero pixels 
in the two maps), divided by the former. The magenta line shows the border of ROI A where the dust fit 
was performed. The maps are shown in the plate carree projection used for the dust and y-ray fits in ROI A. 
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A. 4. Xqo, Dust Specific Opacity and Power in the Local DNM from the 7 - Ray Analysis 

The results in Section can constrain some properties of the local ISM seen in the foreground of the 
HVCs and IVCs, in particular the Xqo = X{}l2)/Wco ratio in the Ursa Major molecular clouds, and the 
average specific power radiated by dust per H atom and unit solid angle, R/A(H), and the average specific 
opacity per H atom, rss 3 /N(R) in the DNM for the other regions. 

These quantities can be extracted from the parameters in Eq. under the assumption that the same CR 
fluxes illuminate the atomic gas traced by Hi and the molecular gas traced by CO or the gas in the DNM. 
Then, Xqo = ^ co / 2 xhi , and R/A(H) = ^ hi/^dnm (when the DNM is traced using R), or r 353 /A(H) = 
^Hi/^DNM (when the DNM is traced using r 353 ). In ROI B we do not detect any significant 7 -ray emission 
associated with the sparsely populated DNM templates, therefore we exclude it from the following discussion. 

First, we check the assumption that the spectrum of CRs illuminating the different phases in the local 
gas complexes is the same in each ROI. Figure shows the xqo and xdnm coefficients as a function of 
energy from the analysis in energy bins in Section |5.3[ Within the large statistical uncertainties the spectra 
of the emissivities in each ROI are consistent with those of Hi in Figure JUL all of them being consistent 
with the spectrum of gas in the local interstellar medium from Casandjian (2014). There is no clear analog 


of the softening in Figure [TT] for ROI A for either CO-traced gas or the DNM. However, the softening may 
be present and hidden by the large statistical uncertainties. 

Given the fact that the spectra are consistent within the statistical uncertainties, we extract from the fit 
coefficients the properties of the local ISM discussed above. We take into account the systematic uncertain¬ 
ties related to inputs to the interstellar emission model using the results of the analysis in Section |5.4.1| and 
from the jackknife tests presented in Section [5. 4. 2[ The ISM properties are reporte d in Table They are 
similar to those inferred from LAT data for other nearby interstellar complexes (e.g., Planck and Fermi Col- 
laborations [2014 ), demonstrating that the modeling of the foregrounds in our analysis is robust. The Xqo 
ratio in the Ursa Major molecular clouds that we obtain is consistent with other estimates from alternative 
methods in the literature (e.g., de Vries et al.||1987 ). 


B. EVALUATION OF ERRORS AND UPPER LIMITS WITH THE JACKKNIFE METHOD 

In this appendix we summarize the method used to evaluate uncertainties and upper limits from the 
jackknife tests in Section 5.4.2 Following Dudewicz & Mishra (1988), given the jackknife estimates 0^ with 

(Bl) 


^ = 1 ,..., let 


1 ^ 
N ^ 


1=1 


be their arithmetic mean. A quantity known as Jj is defined as 

j^ = Ne-{N - i)e,. 

If we define 

1 ^ 

" V F 

the Icr uncertainty of 0 is 


1=1 


N 


The (1 — Of) c.l. upper limit on 0 is 


^ 7V(7V - 1) E 


UL — J + 


(B2) 

(B3) 

(B4) 

(B5) 


where to, jv-i is the upper a*" quantile of the Student’s t-distribution with — 1 degrees of freedom. 
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Fig. 22.— Best-fit 
in ROI A and C. 
uncertainties. The 


Sec. 5.2). 


values for xqo and xdnm (Eq. from the analysis in multiple energy bins (see Sec. |5.3[ ) 
In each panel the points show the values as a function of energy with their statistical 
shaded bands show the Tier interval from the analysis over the entire energy range (see 
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Table 12 

Local ISM properties from the 7-ray analysis 



ROI A 

ROI C 


1.13 ± O.O 8 I 0 73 ± 0.10 


R/Ar(H)’" 

4.1 ± 0.3l° 0 ± 0.3 

3.2 ± 0.3l°;J ± 0.3 

r353/iV(H) = 

2.06 ± 0.06t° “ 



Note. —No CO emission is detected in ROI C. The 7 -ray scal¬ 
ing coefficients for the DNM maps derived from TS 53 in ROI C 
were always consistent with 0 . 

Note. —The statistical uncertainties take into account the co- 
variances of the fit parameters. 

^ 10 ^° cm -2 (K km s“^)“^ 

bio-32 sr-^ 

cm^ R-i 
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